/*
 
 Trabalho 1
 Programação Concorrente
 
 Cálculo de PI com precisão de 10 milhões de dígitos
 Método: Gauss-Legendre (Sequencial)
 Biblioteca de precisão: The GNU Multiple Precision Arithmetic Library (GMP)

 */

#include <stdio.h>
#include <gmp.h>

#define LOG_2_10 3.32192809488736234787

int main (void)
{
    /* Contador de iterações */
    int i = 1;
    /* Variáveis big float da GMP */
    mpf_t pi_old, pi_new;
    mpf_t a_old, a_new;
    mpf_t b_old, b_new;
    mpf_t t_old, t_new;
    mpf_t p_old, p_new;
    mpf_t epsilon;
    /* variável para gravar o nome do arquivo de saída */
    char filename[10];
    /* ponteiro para arquivo de saída */
    FILE * F;
    
    /*
     
     Define a precisão padrão usada nos cálculos:
     
     OBS: Utilizamos 10000100 casas de precisão para
     possíveis overflows durante o cálculo.
     Como o argumento é em bits, e não em casas decimais,
     multiplicamos por uma constante log2(10),
     que é a quantidade média de bits por casa decimal
     
     */
    mpf_set_default_prec(LOG_2_10 * 10000100);
    
    /*
     
     Inicializa as variáveis de acordo com o método:
     a = 1
     b = 1/sqrt(2)
     t = 1/4
     p = 1
     
     */
    mpf_init(pi_old);
    mpf_init(a_old);
    mpf_init(b_old);
    mpf_init(t_old);
    mpf_init(p_old);
    
    mpf_init(pi_new);
    mpf_init_set_d(a_new, 1.0);
    mpf_init_set_d(b_new, 0.5); mpf_sqrt(b_new, b_new);
    mpf_init_set_d(t_new, 0.25);
    mpf_init_set_d(p_new, 1.0);
    
    mpf_init_set_str(epsilon, "1e-10000000", 10);
    
    /* Calcula primeira iteração */
    mpf_add(pi_new, a_new, b_new);
    mpf_pow_ui(pi_new, pi_new, 2);
    mpf_div_ui(pi_new, pi_new, 4);
    mpf_div(pi_new, pi_new, t_new);
    
    /*
     
     Calcula as demais iterações, gravando cada resultado parcial
     em um arquivo de nome it-x.txt
     
     */
    do
    {
        /*
         
         Executa um swap nas variáveis.
         Como a biblioteca trabalha internamente com ponteiros,
         a função é mais eficiente que copiar as variáveis.
         
         */
        mpf_swap(pi_old, pi_new);
        mpf_swap(a_old, a_new);
        mpf_swap(b_old, b_new);
        mpf_swap(t_old, t_new);
        mpf_swap(p_old, p_new);
        
        mpf_add(a_new, a_old, b_old);
        mpf_div_ui(a_new, a_new, 2);
        mpf_mul(b_new, a_old, b_old);
        mpf_sqrt(b_new, b_new);
        mpf_mul_ui(p_new, p_old, 2);
        mpf_sub(t_new, a_old, a_new);
        mpf_pow_ui(t_new, t_new, 2);
        mpf_mul(t_new, t_new, p_old);
        mpf_sub(t_new, t_old, t_new);
        mpf_add(pi_new, a_new, b_new);
        mpf_pow_ui(pi_new, pi_new, 2);
        mpf_div_ui(pi_new, pi_new, 4);
        mpf_div(pi_new, pi_new, t_new);
        mpf_sub(pi_old, pi_new, pi_old);
        mpf_abs(pi_old, pi_old);
        
        sprintf(filename, "it-%d.txt", i++);
        F = fopen(filename, "w");
        gmp_fprintf(F, "%.10000000Ff\n", pi_new);
        fclose(F);
    } while ( mpf_cmp(epsilon, pi_old) < 0 );
    
    /* Libera o espaço de memória alocado para as estruturas */
    mpf_clear(pi_new);
    mpf_clear(pi_old);
    mpf_clear(a_new);
    mpf_clear(a_old);
    mpf_clear(b_new);
    mpf_clear(b_old);
    mpf_clear(t_new);
    mpf_clear(t_old);
    mpf_clear(p_new);
    mpf_clear(p_old);
    mpf_clear(epsilon);
    
    return 0;
}
